Loop Representations of the Quark Determinant 

in Lattice QCD 

A. Duncan^ E. Eichten^, R. Roskies^, and H. Thacker^ 

'^^ • ^Dept. of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15620 

a^ ■ 2pg^j^jlg^b^ p O. Box 500, Batavia, IL 60510 

r^ . ^Dept.of Physics, University of Virginia, Charlottesville, VA 22901 

D 
tin 



> 

in 



c^ 



Abstract 



The modelling of the ultraviolet contributions to the quark determinant in 

lattice QCD in terms of a small number of Wilson loops is examined. Com- 

O I plete Dirac spectra are obtained for sizeable ensembles of SU(3) gauge fields at 

^i^ ' /3 = 5.7 on 6^, 8^ and 10^ lattices allowing for the first time a detailed study 

Q^ . of the volume dependence of the effective loop action generating the quark de- 

^» I terminant. The connection to the hopping parameter expansion is examined 

in the heavy quark limit. We compare the efficiency and accuracy of various 

methods - specifically, Lanczos versus stochastic approaches- for extracting the 

Oh! quark determinant on an ensemble of configurations. 
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1 Introduction 

In a recent paper [|I| we introduced a method for performing an unquenched Monte 
Carlo simulation in lattice QCD in which the infrared and ultraviolet modes of the 
quark fields are treated separately. The low eigenvalues (typically up to a cutoff 
somewhat above Aqcd) are exactly and explicitly calculated and included as a trun- 
cated quark determinant in the the update Boltzmann measure The remaining UV 
modes are included approximately by modelling the high end of the spectrum with 
an effective loop action involving only small Wilson loops. It was found that on small 
lattices (specifically, 6^ lattices at /5=5.7) the accuracy of such a loop fit to the high 
end of the determinant was remarkably good- the x^ per degree of freedom of a fit 
with loops of up to 6 links to the logarithm of the quark determinant, including all 
modes above 340 MeV was 0.23, while the typical excursion of the log determinant 
between uncorrelated configurations was on the order of 20. Of course, one expects 
that the accuracy of such a fit will decrease with increasing volume, and it is not 
clear that this approach will remain practical once lattices of physically useful size 
are reached. As an illustration, we note that recent studies of electromagnetic effects 
in lattice QCD found that the finite volume effects from even the long-range elec- 
tromagnetic effects were controllable on 12^x24 lattices at /5=5.9. A 10^x20 lattice 
at /3=5.7 has almost twice the physical volume, while the lattice discretization effects 
can presumably be substantially reduced by using clover improvement. Our aim in 
this paper is therefore to study the accuracy of effective loop action representations 
to infrared truncated quark determinants for various size lattices at /?=5.7. We shall 
show that reasonably accurate representations of the UV contribution to the quark 
determinant in terms of a small number of Wilson loops are indeed possible on such 
physically useful lattices. To the extent that a small residual error in the loop repre- 
sentation of the ultraviolet part of the quark determinant contributes primarily to an 
overall rescaling of dimensional quantities such as ground-state hadron masses (which 



are dominated by quarks with a limited range of offshellness) such a representation 
should be perfectly adequate for dynamical spectrum calculations on moderately sized 
lattices. 

Past studies of loop representations of the quark determinant P, ^ (which typi- 
cally focussed on the complete determinant, less accurately described by small loops) 
have been hampered by the difficulty of obtaining exact values for the quark determi- 
nant for a sufficiently large sample of independent gauge configurations. Stochastic 
methods |^, |^ can be applied to fairly large lattices but with limited accuracy, whereas 
the direct Lanczos approach loses steam for lattices larger than about 12^. In this 
paper we have employed the Lanczos approach to obtain complete exact Wilson-Dirac 
spectra for sizeable ensembles of 6^, 8^ and 10^ lattices, allowing us to study in detail 
the volume dependence (Section 2) and quark mass dependence (Section 3) of the ef- 
fective loop action fit to the truncated quark determinant at various infrared cutoffs. 
Technical details of the calculations, as well as a comparison of the computational 
burden of the exact Lanczos and stochastic approaches, are presented in Section 4. 

2 Loop Actions for light quarks- volume depen- 
dence 

The main difficulty we encounter in determining the accuracy of a loop representation 
for full or truncated quark determinants in lattice QCD lies in the computational effort 
required to extract complete spectra of the Wilson-Dirac operator for a sufficiently 
large sample of independent gauge configurations on lattices large enough to yield 
physically useful information. A variety of numerical tools, both exact and statistical, 
now exist for accomplishing this task. Technical details of the implementaion of these 
methods will be deferred to Section 4. In this section we present a detailed study 
of the volume and IR-cutoff dependence of loop fits to the quark determinant for 
ensembles of 75 6^ lattices, 75 8^ lattices and 30 10^ lattices at j3 =5.7 and k =0.1685. 



For all of these configurations we have carried out a complete spectral resolution using 
the Lanczos methods discussed in Section 4, and identifying for each configuration 
all 15552 (resp. 49152, 120000) eigenvalues in the 6^ (resp. 8^ 10^) cases. Most of 
these calculations were performed on a 9 node Beowulf system running Linux. 

The configurations used in this study were generated using the truncated determi- 
nant algorithm of [Q. In particular, the update measure included the contributions to 
the quark determinant from low eigenmodes of the hermitian Wilson-Dirac operator 
•y^i^Dw — m) up to a cutoff Aqcd (specifically, we chose a cutoff of 0.45 in lattice 
units, corresponding to about 490 MeV in physical units). This cutoff corresponds to 
the lowest 30 (15 positive and 15 negative) eigenvalues for the 6^ lattices, and to 120 
(resp. 350) low eigenvalues for the 8"^ (resp. 10^) lattices. It is reasonable to expect 
that, by including this infrared contribution in the simulation measure the low energy 
chiral structure is properly treated [Q]. The essence of the task being addressed in 
this paper is then to determine the extent to which the remaining omitted ultraviolet 
modes can accurately be fit by a gauge-invariant loop expansion involving relatively 
few and simple loops. The loop actions discussed here will include only loops (or 
Polyakov lines) with up to 6 links, although the inclusion of loops with 8 links is com- 
pletely straightforward in principle and would of course give a substantial increase in 
the accuracy of the fit. 

On a 6^ lattice with periodic boundary conditions, there are five independent 
gauge invariant contributions to V=lndet{'~f^{Dy[r — m)) involving 6 links or less, 
corresponding to the plaquette (1x1 loop), here denoted Li{U), the length 6 Polyakov 
line traversing the full length of the lattice in any direction, L2{U), and 3 independent 
length 6 closed curves (denoted L3 4 5(^7)), illustrated in Fig.l. In the case of the 8^ 
and 10^ lattices there are only four such terms- the plaquette and the length 6 simply 
connected loops displayed in Fig. 1. In addition the fits contain a constant term 
which can be regarded as the lattice equivalent of the unit operator. We shall be 
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Figure 1: 6 link closed loop contributions to T) 
studying the ultraviolet contribution to T> in which the n\ lowest modes are omitted: 

PK)=Eln(|A,|) (1) 

i>nx 

where Aj are the gauge-invariant eigenvalues of (^^^(Dw — rn)) enumerated in order of 
increasing absolute value. Denoting the approximate loop value of V by Sa and the 
true value by St, the variance per degree of freedom of the fit, a^, is defined as: 



<J^ = Yl i^ain) - St(n)f/{nc - np). 



(2) 



n=l 



for nc configurations and np loop variables (including the constant term). (We note 
here that the terminology "variance per degree of freedom" replaces the usual "chi- 
squared per degree of freedom" as we are dealing with a dimensionless quantity with- 
out statistical errors.) In Figures (2-4) we show the accuracy of the best fit to V{nx), 
for the cutoff corresponding to the actual simulation values - namely nx =30 (resp. 
120, 350 for the 6^ (resp. 8^, 10^) lattices. The a^ is 0.24, 0.28 and 0.79 for the 3 cases 
studied . The very close matching of the loop action to the determinant values for 
0"^ <1 suggests that accurate dynamical calculations should be possible by replacing 
the UV part of the quark determinant by such a loop Ansatz. 

The accuracy of the loop fit to a truncated determinant is increased either by 
including longer loops or by raising the IR cutoff (which requires that more low 
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Figure 2: 6 link fit to ^^so, 6^^ lattices 
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Figure 3: 6 link fit to V^o, 8^ lattices 
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Figure 4: 6 link fit to Psso, 10^ lattices 

eigenmodes be computed explicitly and included in the simulation). The variation 
with cutoff for the three different lattices studied is illustrated in Figure 5. The 
eigenvalue cutoff has been reexpressed as a physical lattice momentum, with the 
conversion performed using an average spectral density obtained by averaging the 
individual spectra for all lattices in the ensembles used. Loops up to length 6, together 
with Polyakov lines stretching across the lattice, have been included in the fit of the 
determinant. The generic behavior is clear from Fig(5). The accuracy of the fit 
improves rapidly with increasing cutoff (initially, the variance decreases roughly like 
Q-Cpmin^ where Pmm = l-^n^l is the eigenvalue cutoff), reaching a small nonzero value. 
For the lattices studied,the minimum attained a^ was less than one in all cases, leading 
to the close matching of actual and loop fit determinant values visible in Figures 2-4. 
The 0"^ then fluctuates more or less randomly for further increases of cutoff around 
this value. This small nonzero contribution at fixed (3 (increasing with the lattice 
volume) is due to higher dimension operators not included in the limited size loops 
in the fit and would presumably be reduced as the lattice /? is increased to push 
the system towards the continuum. The fluctuations are simply a reflection of the 
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Figure 5: a versus Pmin(lat) for 6 ,8 and 10 lattices 

limited statistics of our relatively small ensembles- they become more visible at larger 
cutoff where the signal to background is reduced. This can also be seen by plotting 
the cutoff dependence of the a^ for the 6 link loop fits for the ensemble of 75 8"^ 
lattices broken into three subensembles of 25 configurations each (Fig 6). At smaller 
cutoffs, strong infrared correlations extending over several lattices lead to fairly large 
differences in the a^ of the fits for small Pmin, while at larger values of the cutoff, the 
size of fluctuations in each subensemble is comparable to the difference between the 
cr^ for the different subensembles, suggesting that these fluctuations are statistical in 
origin. 

The accuracy of a simple loop representation for the ultraviolet part of the quark 
determinant, involving the contribution of a relatively small number of Wilson loop 
operators, depends on the fact that only a few independent gauge-invariant operators 
exist of mass dimension 4 or 6 . On the other hand, the number of topologically 
distinguishable Wilson loops grows much more rapidly (in fact, exponentially) with 
the length of the loops. This results in very strong correlations between the values 
of distinct Wilson loop shapes over ensembles of independent lattices. For example. 




Figure 6: a versus p.mm(lat) for 3 ensembles of 25 8 lattices 



the single plaquette value is strongly correlated with the value of the 2x1 loop (L3 
in Figure 1). Thus, although the minimum o"^ attainable for any given cutoff in link 
number is clearly only obtained by employing all Wilson loops up to the prescribed 
length, the actual coefficients of the individual loop values Li exhibit large variations 
from one subensemble to another for a given lattice volume, and also as the lattice 
volume is increased. On the other hand, if only truly independent operators are in- 
cluded, we should expect the coefficents (properly normalized for the overall lattice 
volume) to remain relatively constant as the lattice volume is increased at fixed /3, 
reflecting the contribution of a definite set of low-dimension operators with expecta- 
tion values approaching well-defined values in the infinite volume limit. This can be 
seen in Table 1, where we show the loop coefficients as a function of lattice volume 
for just the single plaquette operator obtained by minimizing the o"^ with respect to a 
fit containing just a constant term (the unit operator) and the single plaquette loops 
(in operator terms, -F^,^). The 10^ lattice results were not included because of the 
limited statistics available in this case (only 30 lattices as compared to 75 lattices for 
the 6^ and 8^ cases). 
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Cutoff (f/a) 
Aqcd 


Volume A^^ 


coefficients 
constant plaquette 


0.00 


6 

8 


-0.0225 
-0.0209 


0.0772 
0.0794 


0.45 


6 

8 


-O.Offf 
-0.0ff8 


0.0695 
0.0699 


0.56 


6 
8 


-0.0043 
-0.0048 


0.0653 
0.0653 



Tabfe 1: Variation with voiume of fit coefficients of fink expansion of tr fn (Det) for various 
fow eigenvafue cutoffs A. 



3 Quark determinant for heavy quarks- hopping 
expansion vs effective loop actions 

The main physical difference between the behavior of the quark; determinants for light 
and heavy quarks lies in the relative variance of the infrared and ultraviolet contri- 
butions. For light quarks, there are important ffuctuations introduced both by the 
infrared modes, which incorporate the proper chiral behavior of the unquenched the- 
ory, and by the ultraviolet modes above Aqcd, which primarily renormalize the scale 
of the theory. Indeed, the latter are quantitatively dominant, but closely matched by 
an effective action involving only plaquettes or 6-link loops which for long distance 
physics reduces to an effective shift in the beta of the simulation. For heavy quarks, 
the density of the quark Dirac spectrum is much reduced in the infrared (which is 
cut off by the large bare quark mass), as are the ffuctuations from the infrared modes 
below Aqcd, while the UV modes still have substantial variance, but again of a form 
which, as we shall show, can be accurately modeled by a simple loop action. The 
relative variance of the low and high end contributions of the (log) quark determinant 
is shown in Fig 7, where the IR/UV cut is placed at the 15 th (positive or negative) 
mode, roughly at Aqcd for a 6'^ lattice at (3=5.7 (for the sake of visibility, an irrele- 
vant constant vertical offset has been applied to bring the various contributions close 
to zero). For the heavier quark, k =0.1500, there is comparatively little variance in 
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Figure 7: IR and UV contributions to quark determinants for light vs heavy quarks 

the IR contributions, but large fluctuations in the UR part. 

The simplest approach one might take for the full quark determinant for heavy 
quarks employs the well-studied hopping parameter expansion for V=\n {■y^^Dw — fn)), 
valid in the limit of small k, i.e. large quark mass. In this section we examine the 
extent to which a truncated hopping parameter expansion can compete with the non- 
perturbative fitting of an effective loop action of the kind described in the preceding 
section. On a 6^ lattice with periodic boundary conditions, we saw previously that 
there are five independent gauge invariant contributions to V of order k® or less, de- 
noted Li,l < i < 5 above. The loop averages Li{U) on a given gauge configuration 
{f/} are normalized to give exactly one on the ordered configuration where all links 
are unity. Then a straightforward combinatoric exercise gives {V= lattice volume) 



V{U) = V{288K'^Li{U)-512K^L2iU)+2304:K^L3{U) 



(3) 



(Note that this computation gives an approximant to the full determinant, with no 
obvious way of implementing an IR/UV cutoff within the hopping parameter expan- 
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Figure 8: Full determinant, hopping parameter vs nonperturbative fit 



sion.) 



Unfortunately, the coefficients of higher order terms in the hopping parameter 
expansion grow fairly rapidly (the number of closed loops increases exponentially 
with the length of the loop) and the expansion converges slowly. For example, for 
fi;=0.15, (with k,^=2.6x10^^ !) the average V over an ensemble of 30 configurations 
gives 172.4, while the hopping expansion gives 147.8. The discrepancy is nevertheless 
mainly due to a few low dimension operators which dominate the contribution of the 
longer loops, as is apparent from Fig 8. The hopping parameter expansion through 
6th order tracks roughly the exact determinants apart from an offset (the identity 
operator). Of course, the nonperturbative fit including Li_5(f/) does much better (the 
a^ is 0.112). We can improve the agreement of the truncated hopping expansion with 
the data by varying both the offset and the single plaquette component from the value 
of the K^ coefficient given in Eq (2). This amounts to including at least the lowest 
nontrivial dimension operator (dimension 4) arising from the longer loops (length 8 
and higher). This ffi is shown in Fig 9, in comparison with the unmodified hopping 
expansion result shifted only by a constant offset. The ffi is certainly improved 
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Figure 9: Hopping expansion fits, with offset vs with offset+plaquette 

by optimizing the single plaquette component (the a^ is now 0.75) but clearly the 
nonperturbative fit of Fig. 8, in which all closed loops through length 6 are optimized, 
still wins by a substantial factor. We can conclude that even for quite heavy quarks, 
the hopping parameter expansion, though analytically available, is not competitive 
with a nonperturbative fit using even a small number of Wilson loops. Determining 
the coefficients in such a fit only requires the extraction of the determinant for a 
few typical configurations. We find that a sufficiently accurate determination of the 
loop coefficients can be obtained by calculating the spectrum on as few as 10 gauge 
configurations. 

We pointed out earlier that the IR portion of the Dirac spectrum is relatively inert 
in the case of heavy quarks. One therefore expects that a nonperturbative fit with a 
few small loops should be accurate for heavy quarks, even if one insists on fitting the 
full determinant, as discussed above. Indeed, we saw above that the a2 with a fit to 
the full determinant is 0.112 including loops up to length 6. The fit is still improved 
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by excising the IR part of the spectrum, though not as dramatically as in the case 
of light quarks where the IR fluctuations require the inclusion of large loops. For 
K=0.15, one finds for example that the cr^ of a 5 loop fit decreases to 0.035, 0.026 and 
0.019 when the lowest 30, 60 or 90 modes are excluded from the determinant. 

4 Extracting the complete Dirac spectrum- explicit 
versus stochastic approaches 

A fairly efficient method for performing a complete spectral resolution of the hermitian 
Wilson-Dirac operator was described some time ago by Kalkreuther |^ . One employs 
the usual Lanczos procedure, without reorthogonalization, pruning out the spurious 



eigenvalues by the Cullum-Willoughby procedure [1^. Typically the extraction of 
the complete Dirac spectrum requires (due to the inexact arithmetic) more Lanczos 
sweeps than the actual dimension (12V, where V is the lattice volume), by a factor of 
2-3. For example, on a 10"^ lattice, the full spectrum is obtained after about 320000 
Lanczos sweeps, as compared to the actual dimension of 120000. The convergence of 
the procedure is improved by starting from gauge-fixed configurations. (In spite of the 
fact that the spectrum is gauge-invariant, it appears that the presence of gauge noise 
can reduce the numerical stability of the Lanczos procedure.) Occasionally, spectral 
fiuctuations lead to two eigenvalues which are almost degenerate, or a real eigenvalue 
almost degenerate with a spurious one, and the spectrum is found to be missing a 
small number of eigenvalues (note that the Lanczos procedure does not identify the 
degeneracy of the various eigenvalues). In an ensemble of 75 6^ lattices, the procedure 
missed a single mode in only two cases. For an ensemble of 20 10^ lattices, 3 eigenval- 
ues were missed 5 times, 2 for 5 configurations, 1 for 6 configurations, and complete 
spectra were obtained for 4 configurations. Although the procedure can probably be 
tuned to reduce the frequency of missed eigenvalues, the general trend is nevertheless 
towards troublesome accidental degeneracies for larger lattices. Recently, we have 
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found that the problem on larger lattices can be considerably ameliorated by a more 
careful retuning of the algorithm for identifying spurious eigenvalues- the resulting 
complete spectra agree completely with the reconstructed eigenvalues discussed be- 
low. Moreover, the convergence of the eigenvalues in the densest part of the spectrum 
occurs very rapidly as one reaches the end of the procedure, so that the spectrum 
of the tridiagonal matrix essentially consists of either spurious eigenvalues (easily 
identified by the Cullum-Willoughby procedure) or accurate converged eigenvalues. 
For example, with 320000 Lanczos sweeps on a 10^ lattice (120000 eigenvalues) one 
finds with a carefully tuned Lanczos calculation exactly 200000 spurious eigenvalues 
and 120000 converged eigenvalues, with the latter satisfying the set of exact sumrules 
discussed below to high accuracy (for a 10^x20 lattice with 240000 eigenvalues, com- 
plete spectra are obtained with 600000 Lanczos sweeps). The present tuning leaves 
about 3 orders of magnitude between the tolerance test for spurious eigenvalues (set 
at 10^^^) and the arithmetic precision (about 10"^"^'^^). For larger lattices where the 
volume and hence the maximum spectral density is a few orders of magnitude higher, 
we expect the Cullum-Willoughby procedure to misidentify some true eigenvalues as 
spurious, leading to incomplete spectra (as we indeed find if the tolerance parameter 
for spurious eigenvalues is set to 10^^°, for example). 

In fact, in all cases mentioned above, incomplete spectra can be completely re- 
paired with the aid of exact sum rules for traces of powers of the Wilson-Dirac matrix 
H = '^^{Dw — m). As we lose at most 3 eigenvalues, the lowest four sum rules suffice 
to determine all missing eigenvalues, with a spare relation left over for checking the 
accuracy of the reconstruction. One easily derives 

(4) 

(5) 

(6) 

Tr(i7^) = 12\/(1 + 64/t2 + (448 - 96 < P >)/t^) (7) 

15 



Tr(iJ) = 


= 


Tr(ij2) = 


= 12V{1 + IQk^ 


Tr(i73) = 


= 



where the Wilson-Dirac operator Dw — m in H is normahzed to be of the form 
1 + ^(hnk terms), and < P > in Eq(7) is the average plaquette value in the given 
configuration. A check of these sum rules on configurations where the entire spectrum 
has been successfully extracted gives Tr(if, H^) <10^^ while the quadratic and quartic 
sum rules are reproduced to at least eight significant figures. Repairing the incomplete 
spectra with these sum rules reveals, as expected, that the missing eigenvalues are 
either close to degenerate with each other or with a converged eigenvalue from the 
Lanczos procedure. 

An alternative to the explicit spectral resolution of H, which is only really feasible 
for small to moderate sized lattices, is the stochastic approach developed by Golub 
and coworkers 0, |[ , and applied by Irving and Sexton in their study [^ of loop actions 
for the quark determinant. Their formalism uses the close connection between the 
Lanczos recursion and Gaussian integration to generate rigorous lower and upper 
bounds to the diagonal matrix element < v\f{A)\v > of any differentiable function 
f{A) of a positive definite matrix A. The spectral sum giving this matrix element is 
transformed to a Riemann-Stieltjes integral and the usual quadrature rules (Gauss, 
Gauss-Radau, Gauss-Lobatto, etc.) applied to this integral can then be reexpressed 
in terms of a Lanczos recursion (for further details, we refer the reader to the paper 
of Bai, Fahey and Golub 0). The use of alternative quadrature rules (in which 
information about the upper and lower limits of the spectrum is included in the 
Gaussian measure) does not seem to matter much in the application to the Wilson- 
Dirac matrix, although we have found that the Gauss-Lobatto version requires about 
50% more Lanczos sweeps to achieve the same precision as the other quadrature rules. 
It is straightforward to generalize the arguments of p to show that the Lanczos 
estimates converge to the correct answer even for non-positive definite hermitian 
matrices (such as the hermitian Wilson-Dirac operator H), although the strict upper 
and lower bounds provided by the formalism no longer hold in the case f{A) = \n{A), 
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as they depend on positivity of the derivatives of f{A) over the entire spectrum. The 
rate of convergence of this procedure is impressive- for an 8^ lattice the estimate of 
< v\ In \H\ \v > for a generic random vector v with only 300 Lanczos sweeps is accurate 
to 7 places on a 8^ lattice, and to about 5 places on a 12'^x24 lattice. 

The Gaussian/Lanczos approach described above immediately leads to a stochas- 
tic method for estimating lndet(|-ff|) = J2i < Vi\ In |-f/^||fj >, where the sum extends 
over a complete orthonormal basis. Namely, one computes an average over a set of 
random vectors \zi > where each vector has components ±1 chosen at random: 

\ndet{\H\) = E{<z\\n\H\\z>) (8) 

vari< z\ In \H\\z >) = 2^|ln|if|./ (9) 

We have checked that the choice of elements ±1 for the components of the random 
vectors is optimal, in the sense that the variance (8) obtained with this choice can- 
not be further reduced by an alternative choice of random variable. This procedure 
therefore gives an unbiased estimator for the full quark determinant, with errors that 
are purely statistical, decreasing as \/yN with the number A^ of random vectors 
used. Evidently the accuracy achieved for a given amount of computational effort 
is directly determined by the size of the off diagonal matrix elements of In | if |. This 
is of course a very complicated functional of the gauge field for a general configura- 
tion. For the ordered configuration, however, H can be explicitly diagonalized, and 
a straightforward computation yields 

var{< z\\n\H\\z >) = 4V{^Y.^n{A{kf + B^{k)B^{k)f 

k 

- {W\n{A{kr + B,{k)B,{k))n (10) 

^ k 

where A{k) = 1 — 2KYl,^cos{k^),B^{k) = 2/T;sin (A;^), and V is the lattice volume. 
In other words, the variance of the estimates (8) is exactly equal to the variance in 
the logarithm of the free quark lattice offshellness (as A(A;)^ + B^{k)B^j_{k) is just the 
lattice version of k^k^ + rn? in the continuum). For not too heavy quarks (it is easy to 
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see that the variance (9) vanishes with k,) the number multiplying the lattice volume 
V in (9) is of order unity, so the optimal stochastic procedure requires A^ ^ V^ to 
determine lndet(|i/|) to within an additive error of order unity, assuming the free 
quark case as a rough guide. Explicit calculations show that the order of magnitude 
of the variance is indeed as given in (9). For example, the free quark formula (9) gives 
about 8x10^ for the variance on an 8^ lattice (choosing k,=0.12), while an actual run 
with 1000 random vectors on a nontrivial 8^ configuration {j3=5.7, /t=0.1685) gives 
a variance of 9.7x10'^, and a final result for the mean (log) determinant is 1224. 0± 
3.1 (the error is obtained by taking the square root of the variance per data point). 
This should be compared with the exact value obtained from the complete spectral 
resolution carried out by the direct Lanczos approach described in the first part of 
this section, which was 1222.148. On a 12^x24 lattice at [3=5.9, k =0.1597, a typical 
configuration gave a variance of 8.0x10*^, compared with 7.2x10*^ from the free quark 
result (9) (again using /t;=0.12 for the free case). An explicit evaluation with 820 
random vectors gave a final result 8682± 9.9 for the (log) determinant in this case. 

We are now in a position to compare the computational efficiency of the direct and 
stochastic methods described above. The stochastic approach yields an estimate of the 
logarithmic determinant accurate to any fixed preassigned error with a computational 
effort growing like V^/"^, while the full Lanczos spectral resolution, which effectively 
determines the determinant to machine precision (actually, about eight significant 
figures), requires an effort of order V"^ . However, the prefactors in each case render 
the Lanczos approach advantageous for small to moderate (say, 12*^) lattices. For 
example, on an 8*^ lattice, the complete spectral resolution requires on the order 
of 130000 applications of the "dslash" (lattice covariant quark derivative) operator, 
while the stochastic method would require about 10000 random vectors, each with 
200 Lanczos sweeps- i.e a total of 2 million dslash operations- to reduce the error 
on the determinant to order unity (four significant figures). On the other hand. 



as discussed previously, on large lattices the direct Lanczos procedure is likely to 
fail as the machine precision will be inadequate to resolve the increasing number of 
accidental degeneracies due to the high spectral density. In this case the stochastic 
approach may be the only option for estimating the complete quark determinant. 
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